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We consider the nonlinear Schrodinger equation (NLSE) in 1+1 dimension with scalar-scalar self 

interaction {tp* ip)'^'^^ in the presence of the external forcing terms of the form re~^^''^^^^ — Sip. 
We find new exact solutions for this problem and show that the solitary wave momentum is conserved 
in a moving frame where Vk = 2fe. These new exact solutions reduce to the constant phase solutions 
of the unforced problem when r — >■ 0. In particular we study the behavior of solitary wave solutions 
in the presence of these external forces in a variational approximation which allows the position, 
momentum, width and phase of these waves to vary in time. We show that the stationary solutions 
of the variational equations include a solution close to the exact one and we study small oscillations 
around all the stationary solutions. We postulate that the dynamical condition for instability is 
that dp{t)/dq(t) < 0, where p{t) is the normalized canonical momentum p{t) — xfp^§^, and q(t) 
is the solitary wave velocity. Here M{t) = J dxtp* {x,t)'il>{x,t). Stability is also studied using a 
"phase portrait" of the soliton, where its dynamics is represented by two-dimensional projections 
of its trajectory in the four-dimensional space of collective coordinates. The criterion for stability 
of a soliton is that its trajectory is a closed single curve with a positive sense of rotation around a 
fixed point. We investigate the accuracy of our variational approximation and these criteria using 
numerical simulations of the NLSE. We find that our criteria work quite well when the magnitude 
of the forcing term is small compared to the amplitude of the unforced solitary wave. In this regime 
the variational approximation captures quite well the behavior of the solitary wave. 

PACS numbers: 05.45.Yv, ll.lO.Lm, 63.20.Pw 



I. INTRODUCTION 



The nonlinear Schrodinger equation (NLSE), with cubic and higher nonlinearity counterparts, is ubiquitous in a 
variety of physical contexts. It has found several applications, including in nonlinear optics where it describes pulse 
propagation in double-doped optical fibers [T] and in Bragg gratings [5], and in Bose-Einstein condensates (BECs) 
where it models condensates with two and three-body interactions [3 H] . Higher order nonlinearities are found in 
the context of Bose gases with hard core interactions [5j and low-dimensional BECs in which quintic nonlinearities 
model three-body interactions [Sj. In nonlinear optics a cubic-quintic NLSE is used as a model for photonic crystals 
[7|. Therefore, it is important to ask the question how will the behavior of these systems change if a forcing term is 
also included in the NLSE. 

The forced nonlinear Schrodinger equation (FNLSE) for an interaction of the form {ip^ip)"^ has been recently studied 
[SI [H] using collective coordinate (CC) methods such as time-dependent variational methods and the generalized 
traveling wave method (GTWM) lO]. In [5] and approximate stationary solutions to the variational solution were 
found and a criterion for the stability of these solutions under small perturbations was developed and compared to 
numerical simulations of the FNLSE. Here we will generalize our previous study to arbitrary nonlinearity {■ip*ip)'^~^^ , 
with a special emphasis on the case k = 1/2. That is, the form of the FNLSE we will consider is 

+ + gir^r^ + S^P^ re-'(k-+'l (1.1) 
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The parameter r corresponds to a plane wave driving term. The parameter 5 arises in discrete versions of the NLSE 
used to model discrete solitons in optical wave guide arrays and is a cavity detuning parameter . We will find that 
having (5 < allows for constant phase solutions of the CC equations. The externally driven NLSE arises in many 
physical situations such as charge density waves [T^, long Josephson junctions T^, optical fibers TJ] and plasmas 
driven by rf fields [15 . What we would like to demonstrate here is that the stability criterion for the FNLSE solitons 
found for k = 1 works for arbitrary k < 2, and that the collective coordinate method works well in predicting the 
behavior of the solitary waves when the forcing parameter r is small compared to amplitude of the unforced solitary 
wave. 

The paper is organized as follows. In Sec. II we show that in a comoving frame where y = x + 2kt, the total 
momentum of the solitary wave Py as well as the energy of the solitary wave is conserved. In Sec. Ill we review 
the exact solitary waves for r — and k arbitrary. We show using Derrick's theorem [TH] that these solutions are 
unstable for k > 2 and arbitrary 6, which we later verify in our numerical simulations. In Sec. IV we find exact 
solutions to the forced problem for r ^ and find both plane wave as well as solitary wave and periodic solutions for 
arbitrary k. We focus mostly on the case k = 1/2. We find both finite energy density as well as finite energy solutions. 
In section V we discuss the collective coordinate approach in the Laboratory frame. We will use the form of the exact 
solution for the unforced problem with time dependent coefficients as the variational ansatz for the traveling wave for 
K < 2. This is a particular example of a collective coordinate (CC) approach [8l flTHlQ] . We will assume the forcing 
term is of the form _ §ip_ t^[\\ choose four collective coordinates, the width parameter /3(t), the position 

q{t) and momentum p{t), as well as the phase of the solitary wave. These CCs are related by the conservation 
of momentum in the comoving frame. We will derive the effective Lagrangian for the collective coordinates and 
determine the equations of motion for arbitrary nonlinearity parameter k. In Sec. VI we show that the equations of 
motion for the collective coordinates simplify in the comoving frame. For arbitrary k we determine the equations of 
motion for the collective coordinates, the stationary solutions {(3 — q — (f> = constant) as well as the linear stability of 
these stationary solutions. We then specialize to the case k = 1/2, and discuss the linear stability of the stationary 
solutions. The real or complex solutions to the small oscillation problem give a local indication of stability of these 
stationary solutions. This analysis will be confirmed in our numerical simulations. For k = 1/2 and r <^ A, where 
A is the amplitude of the solitary wave, we find in general three stationary solutions. Two are near the solutions 
for r = 0, one being stable and another unstable and one is of much smaller amplitude but turns out to be a stable 
solitary wave. 

A more general question of stability for initial conditions having an arbitrary value oi (3(1 — Q) = /3q is provided by 
the phase portrait of the system found by plotting the trajectories of the imaginary vs. the real part of the variational 
wave function starting with an initial value of Pq. These trajectories are closed orbits as shown in [5]. Stability is 
related to whether the orbits show a positive (stable) or negative (unstable) sense of rotation or a mixture (unstable) . 
Another method for discussing stability is to use the dynamical criterion used previously [9j in the study of the case 
K = 1, namely whether the p{q) curve has a branch with negative slope. If this is true, this implies instability. Here 
p{t) is the normalized canonical momentum p(i) = M{t) = J dx4'*{x,t)ip{x,t) and q{t) = v{t) is the velocity 

of the solitary wave. These two approaches (phase portrait and the slope of the p{v) curve) give complementary 
approaches to understanding the behavior of the numerical solutions of the FNLSE. In Sec. VII we discuss how 
damping modifies the equations of motion by including a dissipation function. We find that damping only effects 
the equation of motion for /3 among the CC equations. This damping allows the numerical simulations to find stable 
solitary wave solutions. In Sec. VIII we discuss our methodology for solving the FNLSE. We explain how we extract 
the parameters associated with the collective coordinates from our simulations. We first show that our simulations 
reproduce known results for the unforced problem as well as show that the exact solutions of the forced problem 
we found are metastable. On the other hand the linearly stable stationary solutions to the CC equations are close 
to exact numerical solitary waves of the forced NLSE with time independent widths, only showing small oscillations 
about a constant value of (3. In the numerical solutions of both the PDEs and the CC equations we establish two 
results. First, for small values of the forcing term r the CC equations give an accurate representation of the behavior 
of the width, position and phase of the solitary wave determined by numerically solving the NLSE. Secondly, both 
the phase portrait and p{v) curves allow us to predict the stability or instability of solitary waves that start initially 
as an approximate variational solitary wave of the form of the exact solution to the unforced problem. In Section IX 
we state our main conclusions. 
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II. FORCED NONLINEAR SCHRODINGER EQUATION (FNLSE) 



The action for the FNLSE is given by 



d 



Lett — J dtdx 
As shown in [H] the energy 



E= dx 



K + 1 



is conserved. Varying the action, Eq. (2.1) leads to the equation of motion: 



at ox'^ 



(2.1) 



(2.2) 



(2.3) 



We notice that the equation of motion is invariant under the joint transformation r — >■ — r, ip — >■ — thus if '4'{x, t, r) 
is a solution, then so is —ipix,t,—r) a solution. Letting ip{x,t) = e~^'^^^^^^u{x,t) we obtain the equation 



du(x, t) du(x, t) 
■ 2k 



dt 



dx 



{k^ -6)u{x,t) + 



d'^u{x, t) 
dx"^ 



+ g{u*u)'^u = r. 



Changing variables from x to y where y ^ x + 2kt = x + Vkt we have for u{x, t) — v{y, t) 

.dv^ _ _ ^^^^^^ ^ dMv.t) ^ ^(^.^^.^ ^ ^ 



dt 



(2.4) 



(2.5) 



Note that with our conventions, the mass in the Schrodinger equation obeys 2m = 1, or to = 1/2 so that k = mvk 
Vk/2. This equation in the moving frame can be derived from a related action 



Sy ~ dtdy 



■ *9 



K+1 



(w*w)'*+^ + (fc^ - S)v*v + rv + rv*] 



(2.6) 



Multiplying Eq. (2.5) on the left by v* and adding the complex conjugate, we get an equation for the time evolution 
of the momentum density in the moving frame: py — ^{vv* — v*Vy), namely 



Here 



dpv , di{y,t) ^ 



j{y,t) ^ \ {v^vt - vvt) + \vy\' + {S~ k')\v\' + ^\v\'^^'- 



Integrating over all space we get 



where 



d 
dt 



J dypy{y,t) = F[y = +oo] - F[y = -oo]. 



F[y] = -j{y,t) + rv*{y,t) + rv{y,t). 



(2.7) 



(2., 



(2.9) 



(2.10) 



If the value of the boundary term is the same (or zero) at ?; = ±oo then we find that the momentum Py in the moving 
frame is conserved: 

(2.11) 



Py — constant — j dy-{vv* — v*Vy) 



Using the fact that ■0 = ue and defining 



p= / dx-i^pr.-ri^.), 



(2.12) 
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we obtain the conservation law 

P{t) + M{t)k = Py = constant, 
where M{t) = J dx'ip*ip. If we further define 

Pit) 



P(t) = 



M{t) ' 



then we get the relationship: 



= M{t){p{t) + k)^ constant. 



This equation will be useful when we consider variational approximations for the solution. 
The second action Eq. (2.6) also leads to a conserved energy: 



— / dx v*Vy — {v* v)'"^^ + {k^ — 5)v* V rv ^- rv* 

J [ " K + 1 

Using the connection that ip{x,t) = e^^'^^^^^^v{y,t) we find that E and E^^ are related. Using the fact that 



ipl-ipx — k'^v*v + v*Vy + ik{v*Vy — vVy), 



we obtain 



E^ Ey- 2kPy. 



(2.13) 
(2.14) 
(2.15) 

(2.16) 

(2.17) 
(2.18) 



III. EXACT SOLITARY WAVE SOLUTIONS WHEN 







Before discussing the exact solutions to the forced NLSE, let us review the exact solutions for r = since these will 
be used as our variational trial functions later in the paper. We can obtain the exact solitary wave solutions when 
r = as follows. We let 



ip{x,t) = A sech'^ [I3{x - vt)] e^[p(^-"*)+<^(*)l. 
Demanding we have a solution by matching powers of sech we find: 

p = v/2, 7=1^, A^^^ ^^^^~t^\ ^^^p^ + py^^+s)t + . 



It is useful to connect the amplitude A to {3 and the mass M of the solitary wave. 



42 

M^A^ I dx i;*{f3x)ij{f3x) = —C^; C? 



+00 



dy sech^/''(y) 



(3.1) 



(3.2) 



(3.3) 



One can show that the energy E — dx H of the solitary waves is given by 



E = 



k2(k + 2) 



(3.4) 



For (5 = 0, K > 2 these solitary waves are unstable [HI [19]. For k, — 2 the solitary wave is a critical one. In this 
case, the energy is independent of the width of the solitary wave. Further, in the rest frame w = 0, the energy is zero 
when ^ = 0. When 6 is not zero there is a special constant phase solution of the form Eq. (3.1 1. The condition for 
to be independent of t is 



{p^+p^/n^ + 5) = Q, 

which is only possible for negative 5. For this solution, when p = we find the relationship: 
This solution will be the r = limit of some of the solutions we will find below. 



(3.5) 
(3.6) 
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A. stability when r — 



For the unforced NLSE we can use the scahng argument of Derrick [16] to determine if the solutions are unstable 
to scale transformation. We have discussed this argument in our paper on the nonlinear Dirac equation i33j. The 
Hamiltonian is given by 

H = J dxi^^r.^p. - sr^ - ;^(^»"^'} ■ (3-7) 

It is well known that using stability with respect to scale transformation to understand domains of stability applies 
to this type of Hamiltonian. This Hamiltonian can be written as 

H = Hi-SH2-H3, (3.8) 

where Hi > 0. Here S can have either sign. If we make a scale transformation on the solution which preserves the 
mass M = J ip^t^dx, 

^ p^/^i^(fix) , (3.9) 

we obtain 

H = /3^Hi - SH2 - /3"ff3 . (3.10) 

The first derivative is: dH /d(3 = 2/3Hi~ Kf3'^~^ H3. Setting the derivative to zero at /? = 1 gives an equation consistent 
with the equations of motion: 

Ki73 = 2Hi. (3.11) 
The second derivative at (3 — 1 can now be written as 



h{2-k)H3. (3.12) 



The solution is therefore unstable to scale transformations when k > 2. This result is independent of 5. However once 
one adds forcing terms, it is known from the study of the k = I case |8| that the windows of stability as determined 
by the stability curve p{v) as well as by simulation of the FNLSE increase as 5 is chosen to be more negative. In those 
simulations the two methods agreed to within 1%. 

IV. EXACT SOLUTIONS OF THE FORCED NLSE FOR r / 

For r = one can have time dependent phases in the traveling wave solutions, however when r 7^ one is restricted 
to looking for traveling wave solutions with time independent phases. That is if we consider a solution of the form 

i;{x,t) = e-<'^-+'>^f{y); y^x + 2kt, (4.1) 

then / satisfies 

f"{y)-k'^f + 9{rfrf = r; k'' = - 6. (4.2) 
A. Plane wave Solutions 



First let us consider the plane wave solution of Eq.( 2.3 1 



'(/;(a;,t) = aexp [-i(fcx + 6*)] , (4.3) 



or equivalently the constant solution / = a to the Eq. (4.2). This is a solution provided a satisfies the equation: (here 
a can be positive or negative) 

r = g{a^)''a-{k^ -6)a. (4.4) 



6 



This solution has finite energy density, but not finite energy in general. The energy density H is given by 

g=^(^-+;)f"^^(fe^-.)a^ (4.5) 

Since = 0, the energy is the same in the comoving frame. This is the lowest energy solution for unrestricted g, a, 
fc^ , and 5 The solitary wave solutions we will discuss below will have finite energy with respect to this "ground state" 
energy. There is a special zero energy solution that is important. This solution has the restriction that: 

^^Wf^ =ik^-6). (4.6) 
B. 2k = integer 

When 2k = integer then one also can simply find solutions. The differential equation becomes 

f"-{k'^)f + gir.fr f = r. (4.7) 

Note that again this equation is invariant under the combined transformation / — — / and r — >■ —r. Thus we look 
for solutions of the form f± = ±(o + u{y)) with a > 0; u{y) > so that for this ansatz we have 

\a + u{y)\=a + u{y). (4.8) 

It is sufficient to consider /+ = a + u{y) and generate the second solution by symmetry (r — — r; a — >■ —a; u — >■ —u). 
For /+ we have letting 2k + 1 = A'' = integer 

r = -k''^a + ga^, (4.9) 

and 

N / at\ 

u" = (fc'2 - gNa''-')u - g ^ ^"a^-" . (4.10) 

m=2 ^ ^ 

Integrating once and setting the integration constant to zero to keep the energy of the solitary wave relative to a 
plane wave finite, we obtain 

This equation will have a solution as long as a = k'^ — gNa^~^ > 0. We now discuss solutions to this equation when 
AT = 2, (k = 1/2) and AT = 3, (k = 1). 

C. K = 1/2 

For K = 1/2, N = 2, the equation we need to solve for /_|_ is 

{u'f = au^ -2gau^ - gu^/2; a = k'^-2ga. (4.12) 
This equation has a solution of the form 

u{y) = b sech^^'^iPy), (4.13) 

so that f{y) has a solution of the form: 

f{y)=a + bsech''p{y,t), (4.14) 
where both a,b > 0, and / obeys the equation: 



7 



The wave function is given by 



f" -{k")f + 9\.f\f = r. 



^ exp[-i{kx + 0)][a + 6sech^/3(a; + 2kt)]. , 



(4.15) 
(4.16) 



First consider the case that a > and 6 > 0, so that |/| — f. Substituting the ansatz Eq. (4.14) into Eq. (4.15) we 
obtain 

a^g + 2abgsech^{Py) ~ a [k'f + gsedv' {Py) - b {k'f sedi^i^y) - 6bp'^sech^{/3y) + Ab^'^sech^ (^y) - r = 0. (4.17) 



Equating powers of sech we obtain the conditions: 

- k'^a + ga^ = r; 4/3^ - k'"^ + 25a = 0; - Bfe/J^ + gb^ = 0. 

Or, equivalent ly 

The assumption a > and 6 > requires for consistency that k'^ > and r < 0. So wc can rewrite: 

:i2 



k'^ - + Arg ^2 1 ^— — , 6/3^ 

a= ^ = -Vk'* + Arg, b^^ 



^^^^ , P =^V^' -4|r|5, 
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We also have the restriction that 



2ga > 0. 



(4.18) 



(4.19) 



(4.20) 



(4.21) 



The energy density corresponding to this solution is easily calculated and is given by 
4.ga^ 



H = 



- [k^ - 5)a^ + b^[k^ -6- 2ga + (2/3)6g]sech''[/3(j; + 2kt)] - (4/3)g6^sech''[;9(a; + 2kt)]. (4.22) 



Observe that the constant term is exactly the same as the energy density of the solution (4.3) at k = 1/2. Hence the 
energy of the pulse solution (over and above that of the solution (4.3)) is given by 



E 



384/3^ 
5.g2 



(4.23) 



Because of symmetry there is also a solution with / — > — /, r — > — r. so that we can write the two solutions connected 
by this symmetry as follows: 



f{y) = -sign(r) [a + bsech^ {I3y) ,] 



with a, /? and b given by Eq. (4.20). Thus the wave function can be written as 

'ip{x,t) = -sign(r) [a + bsech^/Six + 2kt)] e-*('=^+^). 



(4.24) 



(4.25) 



Notice that the boundary conditions (BC) on this solution are that for ^^{x) at ±00 the solution goes to the plane 
wave solution: 



ip{x = ±00) — ^ — sign(r)a exp[— i(fca; + 9)]. 
Thus the BC on solving this equation numerically is the mixed boundary condition: 

ikip{x — ±00, t) + 4'xix — ±(X)),t) — 0. 
On the other hand if we go into the y frame where 

u{y — ±co,t) — >■ — sign(r)a, 



(4.26) 
(4.27) 
(4.28) 
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then the BC for the u equation is 



Uy{y — ±00, t) = 0. 



(4.29) 



Consider the case where g — 2,k = —0.1,6 = — l,r = —.01. For this choice of parameters the exact solution is 
given by 



f{y) = 0.727191sech2(0.492338y) + 0.0101031. 



(4.30) 



The stationary solution found using our variational ansatz which we will discuss below (See Eq.(6.66l) yields the 
approximate solution: 

/2+(y) = 0.745042sech2(0.498345y), (4.31) 
which is a reasonable representation of the exact solution as seen in Fig. [T] We will show later by numerical simulation, 




FIG. 1: Exact f{y) (solid line) versus variational solution f2+{y) (dashed) for k — —0.1, 5 — —1, g = 2, r = —0.01. 

that this solution is metastable in that f3 remains constant for a short period of time and then oscillates. In the 
numerical simulation of Fig. [9]the parameters used in the simulation are k = 1/2, 5 = —1, g = 2,k = 0.1, r = —0.075. 
For that case: 



f{y) = 0.486113sech^(0.402539y) + 0.0904622. 



(4.32) 



Note that here r is 14 % of A so that the variational solution is worse. The unstable stationary variational solution 
corresponding to this is 



/2+(2/) = 0.745535sech^(0.498509y). 
To make a comparison we plot the square of these solutions in Fig. [2] 



(4.33) 



D. Finite energy solutions 



FromEq. (4.22 



we notice that we can have a finite energy solution (with energy as given by Eq. (4.23 1) when the 
constant term in the energy density is zero. This leads to the condition Aga = 3fc'^, and yields b = —a. This solution 
is not a small perturbation on the r — solution which has a = 0, 6 ^ 0. Instead for the finite energy solution a — —b 
and the form of the solution is 



/(y)=Atanh2 /3y, 



(4.34) 



which has the appearance of Fi g. [S] 



If we insert the solution Eq. (4.34 1 into Eq. (4.15), we obtain 



Ag\A\ tanh*(/32/) + 2^/3^ - Ak'^ tanh^(/3y) + 6A^^ tanh*(/3y) - 8A/3^ tanh^(/?2/) - r = . 



(4.35) 



9 




FIG. 2; Exact f^{y) (solid line) versus variational solution fi+iv) (dashed line) for 5 — —1, g — 2, k = 0.1, r = —0.075. These 
parameters were used in the simulations shown in Fig. [9] 



Thus Eq. (4.34) is a solution provided 



(4.36) 



which requires that both g < and k'^ < 0. We see this solution has the symmetry that r changes sign with A and 
thus the sign of / depends on the sign of r. Because the coupling constant needs to be negative, the quantum version 
of this theory would not have a stable vacuum. We can write this solution as 



— tanh liy, 
\r\ 191 



(4.37) 



where now = \k'^\/8 and r has the special values: r± — ±{3/16){k'^)'^ /\g\ 




FIG. 3: Finite energy solution for k = 1/2,/3 = 1, normalized to unit height. 



For K = 



1. Periodic Solutions for k = 1/2 

1/2 one can easily verify that there are periodic solutions of Eq. (4.2) of the form 

/+(y) b dn^(/3 y,rn) + a. 



(4.38) 
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where dn{x,m) is the Jacobi Elhptic Function (JEF) with modulus m. Again for a > 0; , 6 > we have |/+| = 
Matching coefficients of powers of dn one obtains 



_ \/l - m + m2fc'2 - (2 - m)A/4gr + k"^ 
2g\/l — m + 

5 ' 4\/l - m + m2 

The requirement that a > translates into a restriction on r, namely r < —3(1 — m)k'^ /4:g{2 — m)^ 



(4.39) 



E. K,= l 



When K — 1, the equation for / is 



r-k''f + g\.f\\f = r. 



(4.40) 



Again it is only necessary to consider = (a + u{y)) with both a > and u{y) > 0. For k = 1, = 3, the equation 
we need to solve for u{y) is 



{u'f = av? - Igav? - gu^ 12. 



This equation has two solutions of the form 



c ± cosh(/32/) ' 



(note y — x + 2kt). These solutions were found earlier by Barashenkov and collaborators, 
u for /+ requires that > 0, c > and we choose the + solution. This leads to 



7/^+ (a;, t) = exp[— i(A;a; + 9)] 



provided 



f3 = Va; b 



V2c 



c + cosh 13 {x + 2kt) 
\/2ag 



ga 



; a = k — 3ga 



(4.41) 
(4.42) 

Our condition on 
(4.43) 

(4.44) 



Since a > 0, we require k''^ > 3ga^ > 0. Here r = —k'^a + ga^ < —2ga^is negative. When we let a — >■ —a, b — > —b, c 
— c, then r changes sign and becomes positive. Thus we can write 



V'i (a;, t) = — sign(r) exp[— i(kx + 0)] 
These solutions are non-singular since 



|b| 



|c| +cosh/3(x + 2kt) 



1 - = ga/[ga + 2g^a'^] > . 
The energy density corresponding to these solutions is given by 



H 



3ga^ ,,2 2 262(fc'2 -3ga2) 29{cp^+gab) 9[I3^{^ - I) - {g/2)b^] 



k''a' + 



[c + cosh(/3?/)]2 [c + cosh(/3y)]3 



[c + cosh(/3y)]4 



(4.45) 



(4.46) 



(4.47) 



Note that the constant term is exactly the same as the energy of the solution (|4. 3[) as given by Eq. (|4.5|) at n = 1. 

;n by 

(4.48) 



Hence the energy of the k, — 1 pulse solution (over and above that of the solution (4.3)) is again finite and given by 

8gay^ 



E = 



8a'^/2 lQ^/2a'^/'^ag 

•-'2 ^ T , „ o 9iQ/9 -'3 



gaA-2a?g'^ [ga + 2a'^g'^Y'/'^ [ga + 2a^g^]^ 
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where 



dy 



We have 



dx 



[c + cosh(/3y)]J 



tan 



- + b cosh(a;) y/b^ - a? 



J =2,3,4. 



- a? 
b + a 



, if b^>a^ 



(4.49) 



(4.50) 



Note that in our case a = c,b ~ l,c^ < 1. From here it is easy to calculate the integrals I2, ^3, ^4 and hence show that 
the energy of the + and the — solution (over and above the solution (4.3)) is given by 



8«i/2 2 8V2a 
— [a + 3a g) — 

35 Vff 



{a + 20^ g) tan ^ 



2a^g 2a^g 



a 



(4.51) 



1. Periodic solutions for k = 1 

For K = 1 we can generalize the solitary wave solution for 

r-k\f + gf = r, 

namely 

b {ac + b)sech[dy] + a 



f = a + 



c + cosh[(iy] 1 + c sech[dy] 



(4.52) 



(4.53) 



and obtain a periodic solution in terms of the Jacobi elliptic functions. The generalization of sech[c??/] is the Jacobi 
elliptic function dn(d?/,m) which also has the property |dn| = dn{dy,m) 
One finds that 



/+ 



a + h dn{dy, m) 
1 + c dn{dy, m) 



with a > 0\ ,h > 0,c > obeys |/+| = /+ and is an exact solution to Eq. (4.52) provided 



{ac + h) [agh — ck'"^) 



2c2 



iagh - ck'^ 
2c{m - 2) 



(4.54) 



(4.55) 



Ad^ 



, (1 — m)h 



2 _ 1 + 2aga^ - k'^a 
Ad'^a'^g'^a?- 



.g[2 + a.ga2(l-fc'2a)] 
Here a = c/gah which obeys a cubic equation: 

16(1 - m)a^d'^ga^ = [2 + aga^{l - k''^a)][l + 2aga^ - k'^a] 



(4.56) 



(4.57) 



F. K = 3/2 

For K = 3/2,A^ = 4we have instead for /+ that u obeys the equation 



{u'Y = au^ - Aga^u^ - 2gau^ - 251*75. 



For this case, the formal solution: 



y + c= dy 



yy/a- Aga^y - 2gay^ - 2gy^/5 



(4.58) 



(4.59) 



leads to an elliptic function. 
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V. VARIATIONAL APPROACH TO THE FORCED NLSE 

For the problem of small perturbations to the unforced problem we are interested in variational trial wave functions 
of the form: 

V.i(x,t) =^(t)/[/3,(i)(x-g(i)]e^tP(*)("-'(*»+^(*)l. (5.1) 
Here we assume that the collective coordinates (CCs) A{t), (l){t),p{t), q{t) are real functions of time and that / is real. 



On substituting this trial function in Eq. ( 2.1 ) and computing various integrals one finds that the effective Lagrangian 
is given by 



L = C,^\ p{t)q{t) - m - {p\t) -5)- -^Pt 

„[ /1('+^12(k+1) 

n - 4r cos[kq{t) + m+ 0] IW) + fc, , (5.2) 

Py{K + 1) 



where 



while 



oo 



Cu = / dy[f{y)Y(-+^^ , D,^ [/'(y)]^ , (5.3) 



I[p{t) + k],(3,] = / dyfiM cos[(p(t) + k)y] . (5.4) 





Note that for our parametrization of the variational ansatz: 

rMt) = i^ird:,ij-^d,r)- (5.5) 

Thus from our previous discussion about the conservation of momentum in the comoving frame with k = v/2 we have 
that p{t) and M(t) are in general not independent but satisfy 

M{t){p{t) + k)= constant. (5.6) 

For the case where p{t) + fc 7^ 0, one then gets the relation 

M{t) = Constant/(p(t) + k). (5.7) 

A. Variational Ansatz 

Generalizing the k = 1 choice of [9] for / to arbitrary k, we choose 

=sechi/'^[y], (5.8) 

so that our variational ansatz is 

V'„i(a;,t) = A(i)sechi/"[^,(i)(x-(7(i)]e'[P(*)("-«(*»+^(*)l, (5.9) 
where Ait) is of the same form as in the unforced solution, namely 



/3i/-yaM; a{n)^(^^^ . (5.10) 



Note that we have chosen / to be real. In our previous studies of blowup in the NLSE we used a more complicated 
variational ansatz where there is another variational parameter in the phase multiplying the quadratic term [x — q{t)Y 
[IS] [TO] . Defining the "mass" of the solitary wave M{t) via 

M{t) = / dxV'*^' = CoA{tf/Py{t) = CoP'^'''-^a{K), (5.11) 
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we then find for the ansatz Eq. (5.10) that for p + A; 7^ 

+ k) = constant. (5.12) 
For the forcing term contribution to the Lagrangian, Eq. (5.2), we need the integral 

noo 

/[i^,a,/?„]= / dxcosiax^sech^i/Syx) ^[2''-^/f3^r{i')]r{iy/2 + ia/2/3y)r{iy/2~ia/2f3^) (5.13) 
Jo 

where Re/3y > 0, Re{v) > 0, a > 0. Special cases of this are obtained for = 1, 2, 3. We find 

TTSCchf^) Tracschf^) 7^ (a^ + fl^) sech 

= /[2,a,/3.] = /[3, a, = " ^° 4^3 (5-14) 



For our variational ansatz Eq. (5.8) 



so that 



y^r (i) v^r (1 + i) v^r (i) 



Thus we can write everything in terms of Cq. 

B. Arbitrary k 

For arbitrary k the effective Lagrangian is 

L M{t) [p{t)m - m ~ {p\t) -5)+ 

2i r /3y''"V"('«)r+r- cos{kq{t) + (t){t) + 6) 

where 



r± = r 

Introducing the notation 



(5.17) 



(K^W^^;))^ «(.)^«'-., MS) 

B = ^ + kq + 9, (5.19) 
an important special solution is obtained in the limit C — )• 0, i.e. p{t) — —k. In this case the Lagrangian becomes 



L[C = 0] = M{t) (^-kq{t) - 0(t) - (fc2 -S)+ 
2i r /3y'^-^y^r2(2^)cos(i?) 



(5.20) 



This leads to the equations: 

q = -2k, (5.21) 



I 

(2--)rM^) 



_ .2 , ^ , 2^ ra;:^/^;3-^/-(i-.)r(i + i)r^(^)cos(i?) ^ 
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Assuming the ansatz /3 = f3s and 4>{t) = <t>s — ctgt in Eqs. ( 5.22 )-( 5.23), where f3s, as and (j)s are constant, we obtain 

as = -2fc^, (f>s = nn - kqo - 9, (5.24) 

where n is an integer and /3s sohition of 

These equations become much simpler in the comoving frame, so we will now turn our attention to solving the 
problem in that frame. 

VI. VARIATIONAL APPROACH FOR THE ACTION IN THE COMOVING FRAME 

In the comoving frame we use the following ansatz for the variational trial wave function: 

= A(0/[^.(t)(y-9(t)]e*[^(*^(^-'^(*»+^(*)]. (6.1) 
Comparing with our general ansatz we have the relations: 

y ^ x + 2kt; p^p + k,q = q + 2kt, and 4) = (f> + kq + 0. (6.2) 



On substituting these relations in the effective Lagrangian ( 5.2 ), the Lagrangian in the comoving frame takes the form 



R \ \ n " cos[0(t)]/[p(^), /?]. (6.3) 

We notice that the canonical momenta to and q are simply given by 

^ = M{t)p{t): ^ = ~M{t). (6.4) 

A. Variational Ansatz function choice 

Again choosing from Eq. ( |5.8[ ) our variational ansatz for v is 

v,aiy,t) = A(i)sechi/'^[/3,(i) {y - g(i)] e'[pW(y-m)+m , (6.5) 



where A{t) is again given by Eq. ( |5.10 1. Again the "mass" of the solitary wave M{t) is given by M(t) — J dxv*v — 



CoAitf / I3y{t) = Co/3--'^a{K) and from the conservation of = Mp (Eq. plsj )) we find for p 7^ 



(3" "'^p(t) = constant. (6.6) 



We can rewrite Eq. (6.3) as 

L = M{t) I p{t)q{t) ~ c^{t) - {f{t) + k'-5)+ 0t 



2/j.\ I ;2 r\ I o2 (^ '*) 



k2(2 + k) 

2^ r /3y''"^y^r+r_cos(^(0) 

where r± = F (^T^ + k)) • Lagrange's equation for q we obtain 

^ [M(t)p(t)] = 0, ^ l3'^/''-^p = P^''^~^[p{t) + k]= constant. 



(6.7) 
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Letting 



so that 



G{x^P/p) = 



2i r ^a(K)r+r_ 



dG_G^ dG_ pG' 



we get the following differential equations: 



and for 6 we have 



.^ = 2p+^G'/3V^-eos^W, 



{2-k) 



^ (p{t)m - 4>{t) - (f(t) + k^- s)) 



2 (2 - ^) /3^/" f P ^, l~>i 
k3 + M{t) \(3^ 



G' —G 1 cos4). 

up 



(6.9) 

(6.10) 

(6.11) 
(6.12) 



(6.13) 



B. p/p 

Introducing the notation C = 'np/213, a special stationary solution is obtained in the limit C — ^ 0, i.e. p{t) = 0. In 
this case the Lagrangian (6.7) becomes 



L[C = 0] 



/^r(i)a./3.(t)2/-i 



V J- 



5 - m ~k'+ K 



2 I o2 (2 '^) 



k2(2 + k) 



2i r ^,pI''-'t\^)cqs{^) 



This leads to the equations: 



g = 0, 



- r A — sin(0), 

A hi 

5-e + r 



1 - K 

2 - K 



C0S((^), 



where 



2t ^a'^'l'T{\ + l)T^{lJ 
^/^^) ' 



A = 



B = 



A > 0, 



(6.14) 

(6.15) 
(6.16) 

(6.17) 



Assuming the ansatz fi = fig and 4){t) — (ps in Eqs. (6.16)-(6.17), where Ps, and (ps are constant, we obtain 



where n is an integer and /3s the solution of 



{2-k) 



= 0. 



(6.18) 



(6.19) 



(6.20) 
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1. Linear Stability 



Let us look at the problem of keeping q = p = Q, and looking at small perturbations about the stationary solutions 
/3 = (is and = = 0,7r. The analysis depends on whether rcost^g > or rcos^s < 0. We discuss the two cases 
separately. We let p = |rcos0s| > 0. 

Case lA: r cos = p > 

The linearized equations for 5/3, 5(j) then become 



(k-1)/k 



-CiC 



(6.21) 



5^ = 



2(3s p B{1- K 



k(2 - k)A 



6(3 = C2S13 . 



(6.22) 



Combining these two equations by taking one more derivative we find: 



< + Sf5(t) = 0, 



(6.23) 



where fi^ = C1C2. Whether this will correspond to a stable solution will depend on signs of c\,C2 and hence ciC2- It 
turns out that the answer depends on the value of k. In particular, the answer depends on whether k < 1 or 1 < k < 2 
or K > 2. So let us discuss all three cases one by one. Note that the above analysis is only valid if k ^ 2. If k = 2, 
the entire analysis needs to be redone. 

K < 1: In this case ci, C2 > and hence $1^ > so that the solution is a stable one. 

1 < k; < 2: In this case, while ci > 0, the sign of C2 will depend on the values of the parameters. In particular, if k 
is sufficiently close to (but greater than) one, then C2 > and hence fi^ > so that one has a stable solution. On the 
other hand if /3 is sufficiently close to (but less than) two, then C2 < and hence fJ^ < so that one has an unstable 
solution. In particular, no matter what the values of the parameters are, as k increases from one to two, the solution 
will change from a stable to an unstable solution. 

2 < k: In this case, while Ci < 0, C2 > and hence < so that one has an unstable solution. 
Case IB: r cos 0s = —p <0. 

The linearized equations for 5p, Scj) then become 



(6.24) 



5^ = 



2^ _ pB{l-K) 



6(3 = . 



(6.25) 



Combining these two equations by taking one more derivative we find: 



6^ - n^6p = 64>- n^64> = o. 



(6.26) 



where fl^ = C1C3. Whether this will correspond to a stable solution will depend on signs of ci,C3 and hence C1C3. 
Again it turns out that the answer depends on whether k < 1 or 1 < k; < 2 or «; > 2. So let us discuss all three cases 
one by one. 

K < 1: In this case Ci > while the sign of C3 and hence fl^ depends on the values of the parameters. For example 
if K is sufficiently close to one, then C3 > so that the solution is an unstable one. 

1 < K < 2: In this case, both Ci, C3 > 0, and hence > so that the solution is an unstable one. 

2 < k: In this case, while ci < 0, the value of C3 will depend on the value of the parameters. In particular if k 
is sufficiently close to two then C3 < and hence $1^ > so that one has a stable solution. On the other hand for 
K > Kc > 2, C3 > and hence < so that one has an unstable solution. In particular, no matter what the values 
of the parameters are, as k increases from two, the solution will change from a stable to an unstable solution. 

We now discuss the k= 1/2 case in some detail. 
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C. K = l/2 

For K = 1/2 the effective Lagrangian is 

P (^-ngrmcsch cos(^(i)) + 4/3(t)3 [p{t)m - Pit? - h) + 5 ~ e) + f • (6.27) 

From the Euler-Lagrange equations we find that: 

I3^{t)p{t) = constant, (6.28) 



and further, 



■ 5rC2 sin 

^^^^^h^hc' (^-^^^ 

7 ~2 , .o2 , e ,2 , .gTrrpcos^ grC^ coth C COS 

+4/3 +5-fc + 4/33 3i,hc [^-^^"^^^^- 6/3^sinhC ' ^^-^'^ 

We are interested in the stationary solutions of these equations. We assume 

q = Vst, P = Ps, P = Ps, 4> = 4>s- (6.33) 

We now have Cg = j^iPs) and we need sin 0s = so that 0^ = 0, vr and cos (jjs ~ ±1. 
The equations for the two choices of cos (ps become: 

V, = 2ps ± , ^/^'^T^ ^ [1 - a coth C] , (6.34) 
4p^ smh Cs 

0=Ps+4/3s-fc ± 4^3 3i,hcJ ^'^-^°^^^-]^ 6/3gsinhC. ' ^^-^'^ 
This leads to two families of stationary solutions on two curves in the three dimensional parameter space of Vs,/3s, 



and Ps- Equation (6.351 can also be written: 



1. phase portrait 

The variational wave function for k = 1/2 in the comoving frame is given by 

v,{y,t) = ^sech^/? [y - q{t)] e'{p^*^^y-'i^*'>+M . (6.37) 
9 

The position of the solitary wave q{t) consists of a linear term vt plus an oscillating term. The average velocity v 



can be obtained from the numerical solution of the ODEs (6.29)-(6.32) for the collective variables. The relevant wave 



function for the phase portrait ^ is to evaluate Vy at the point y — vt . The phase portrait is obtained by plotting 
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the real versus imaginary part of the resulting wave function as a function of time by solving the ODEs for various 
initial values of /?. That is, we consider 

^ = ^sech2/3(0 [vt - q{t)] e4p(*)(«**-9(*))+'^(*)] , (6.38) 
9 

and plot Im ^! vs. Re 5* for fixed parameters varying the value of /3q. Orbits which are ellipses in the positive 
(negative) sense of rotation predict stable (unstable) solitary waves in the simulations. When the orbit has both 
senses of rotation, as in a horseshoe shape, the solitary wave is unstable. These behaviors are shown in the phase 
portrait of Fig. |4] The three stationary solutions that we found earlier correspond to the fixed points of the phase 
portrait. The stability of these fixed points can be determined either numerically by solving the CC equations, or 
analytically by a linear stability analysis that we will now discuss. All the fixed points are on the real axis since 
= 0, TT for the stationary solutions. We have at the fixed points 

^^ = Mie"Ae. (6.39) 
9 

We remark that a stable fixed point only means that the CC solutions in its neighborhood are stable, but an orbit 
around this fixed point docs not necessarily have a positive sense of rotation. An example is the medium size ellipse 
in Fig.|4} 




FIG. 4: Phase portrait; imaginary vs. real part of the wave function Eq. (6.38 1. Orbits with positive sense of rotation predict 
stable solitons in the simulations. If the orbit, or part of it, has a negative sense, the soliton is unstable. The filled and open 
circles are stable and unstable fixed points, respectively, Eq. ( |6.39[ ). The orbits are obtained for Po — 0.2 (medium size ellipse), 
Po = 0.53 (small ellipse), /3o = 0.63 (horseshoe) and /3o = 0.7 (large ellipse), keeping fixed 00 = 0, po = —k + 10~^, qo = 0, 
except for 0o = tt and /?o ~ 0.498345 (separatrix, red curve), see Fig. 18 The parameters are r — 0.05, k — —0.1, S = — 1, 
en = and a = 0. 
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2. Linear stability analysis of stationary solutions 
To study the stability of these stationary solutions we let 



From Eq. (6.30) we obtain: 



We can schematically write Eq. (6.32) as follows: 

^ = p2 + - k'^ + F[p, /?, C] cos 4), 



where 



4/33 sinh C 

This leads to the equation for linear stability 



F{P,I3,C]= , J.,j:^ [l-gcothC] 



6^2 sinh C 



' dF dF dF 
= 2ps5p + 8/3,^/3 ± ( —5p + — <5/3 + —5C 



(6.40) 



(6.41) 



(6.42) 



(6.43) 



(6.44) 



where the derivatives are evaluated at the stationary values PstCstPs Now the conservation of momentum in the 
comoving frame leads to 



P^p — Ci — constant. 



So we have 



(6.45) 



(6.46) 



Putting this together we get: 



where 



C0 



We can write 



9F...a^(9F ^CsdF 



so taking derivatives of Eqs (6.48) and (6.50) and combining we get 



54> ± cpc^54> = 0, (5/3 ± cpc^5j3 = 0. 
Thus, depending on the sign of C/jc^ we will have oscillating or growing (decreasing) solutions. 



(6.47) 



(6.4 



(6.49) 



(6.50) 



(6.51) 



Once we have solved for (5/3, we can go back and solve for 5q. We can write the Eq. (6.29) for q as 



^= 2p + cos(^ i^i(/3,C), 



(6.52) 
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where 



Letting q = Vg + Sq we have 



f dF dF 
Sq -^SP±[-S^ + -SC 

' U/3 dC 



Cn±Sf3 



Thus if we are in a region of stabihty so that 



(5/3 — ep cos(rit + a), 



then we obtain 



5(3. 



5q = Sq{0) + ^^5^ [sm{nt + a) - sin a] 



(6.53) 



(6.54) 
(6.55) 
(6.56) 



D. Dynamical stability using the stability curve p{v) 

In references [5| [S] it was shown for the case k = 1 that the stabihty of the sohtary wave could be inferred from 
the solution of the CC equations by studying the stability curve p{v), obtained from the parametric representation 
p{t), v(t) ~ q, where p is the normalized canonical momentum in the lab frame and v = q is the velocity in the lab 
frame. We can determine these quantities using the relations 



p{t) = p — k; q = q — 2k. 



(6.57) 



A positive slope of p{v) curve is a necessary condition for the stability of the solitary wave. If a branch of the p{v) 
curve has a negative slope, this is a sufficient condition for instability. In our simulations we will show that this 
criterion agrees with the phase portrait analysis. 



E. K=iC = 



For the special case that C = 0, using the identity x cscLt — > 1, the effective Lagrangian Eq. (6.27) becomes: 



L[C = 0] = -2grP{t) cos(0(t)) + 4/3(i)-^ [6 - k^ - (j){t) + 



(6.58) 



This leads to the equations: 



p = 0;g = G, 



2 , Aa2 gr 
6/32 



cos(0). 



(6.59) 
(6.60) 



The stationary solution in the comoving frame is represented by /3 = /3s , q — qo and 4i{t) = (ps, where 4>s is given 

by 

4>s ~ niT, (6.61) 
with n being an integer. It is sufficient if we take n = 0, 1. For n — and r > or equivalently, n = 1, r < we have 



^k"^ + 8g\r\/3 



(6.62) 
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whereas for r < 0, n = or equivalently r > 0, n = 1 there are two sohitions given by 

Vfcl_85M/3^ fc'V(5H) > 8/3 « 2.66667. (6.63) 

Note when C = then ps — 0. We expect that these stationary solutions are close to exact solitary wave solutions to 
the original partial differential equations that may be stable or unstable. 

We are interested in seeing how these solutions compare to the exact solution we found earlier as well as to the 
unperturbed r = exact solution. The unperturbed solution with r = 0,p = is given by 

/o(2/) = 6^sech2[^y], (6.64) 

where = —(5/4. Thus for (5 = — 1, g = 2 the exact unperturbed solution is 

3 



UV) - ^sech^ (I) . (6.65) 



The exact perturbed solution for r — —.01 and fc — —0.1,(5 — —1, is given by Eq. (4.301 i.e. /(y) = 
0.727191sech2(0.492338y) + 0.010103. 

Let us now look at the three stationary variational solutions. Since r < we have for n = the two possibilities 
based on the choice of the ± in the square root. For the positive choice we obtain: 

/+(y) ^ 0.745535sech^ (0.498345?/), (6.66) 

which is very close to /q. We will find from our linear stability analysis that this is unstable. 
On the other hand, for the negative root we obtain 

/_(y) = 0.011965sech^(0.06315322/), (6.67) 

which is far from /g. This solution is stable to small perturbations. For n = 1 we obtain instead 



/3? = ^ + ^^^!±M^ ^ 0.255148, (6.68) 
8 8 

so that 

h{y) = 0.765443sech2(0.50512l2/). (6.69) 

We find that this is stable to small perturbations. 

We are also interested for comparison purposes in using the same parameters we used in the study of the k — \ 
problem [5], i.e. fc = —.1,(5 — — l,.g = 2,r = 0.05, (70 = 0, (/)s = 0. Then for this positive value of r, the amplitude of 
the n = solution is (note y = .t + 2fct) 

/i (y) = 0.804134sech2 [0.51773y] . (6.70) 

The phase of the solution is zero. This is not very far from the value of /q given in Eq. ( 6.65[ ). The comparison is 
shown in Fig. [5j This solution turns out to be stable under small perturbations. If we instead chose fe = +0.1, the 
result for /i(y) would be the same but the solitary wave would move in the opposite direction. 
The two solutions for rt = 1 = tt are 

/2+(y) ^ 0.704252sech^(0.484511y); h-Wi = 0.053248sech^(0.133227y). (6.71) 

/2+(y) is again similar to /o but /2-(y) is clearly not. Our linear stability analysis leads to the conclusion that 
is unstable but f2-{y) is stable. In fact, /2+(y) is unstable to linear perturbations but is actually mctastable and the 
original solution switches to a separatrix solution at late times. /2- is stable to small perturbations. The separatrix 
for these values of the parameters is shown below in Fig. |18| 

Next we want to consider initial conditions when r — 0.005 and k = —0.01,(5 — l,.g = 2. For these initial 
conditions it is easier to deal with the periodic boundary conditions on the numerical solutions of the PDEs. The 
n = solution is 

/i(y) = 0.755042sech^(0.501678y) , (6.72) 
and we find that this is stable to linear perturbations. The positive and negative roots for the n = 1 solution are 

/2+(y) = 0.745042sech^(0.498345y); /2_(y) = 0.00503328sech2(0.0409604y). (6.73) 
We will find using linear stability analysis that /2+(y) is unstable but /2-(y) is stable. 
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FIG. 5: foiu) vs. forced variational solution for fi(y) for k — —0.1, S = —1, g = 2, r — 0.05. 



F. linear stability at k = |, C = 



Let us look at the problem of keeping q — p — 0, and consider at small perturbation of Eqs. (6.59) about the 
stationary solution (3 — (5^ and (j) — (ps — Ojir. When (jig = 0, the linearized equations for Sf3, 5(j) become 

5P = -^5i^-cM. (6.74) 

Here the sign of ci is given by the sign of r. 

H = [^Ps + 5p = C25I3 C2 > 0. (6.75) 

Combining these two equations by taking one more derivative we find: 

5'l3 + Q^SP = + if 54) = 0, (6.76) 

where = ciC2. Thus, for the case when we have an exact solution and we choose r = — 0.01, = 2, fc = —0.1, 5 = —1 
we find that for the stationary solution f+{y) we obtain ci = —0.0066866, C2 = 3.93426, and ft^ = —0.0263068, 
suggesting an unstable solution. For /_(j/) we instead obtain ci = —0.0527817, C2 — —2.56198, and fi^ — 0.135226, 
suggesting a stable solution. When r > 0, the solution fi{y) which corresponds to r > 0, (/> = is always stable. 
If instead we look at the small perturbations around the solution where 0^ = tt, we instead obtain 

= ^S4> = ciS4>. (6.77) 



Here again the sign of ci is determined by the sign of r. 

U = (813, - 1^3^ 5/3 = 0^5/3. (6.78) 

When r < and we choose the previous values- r = —0.01, g — 2,k — —0.1, S = —1, k'^ ~ 1.01 we obtain for fi{y) 

Pi ^ 0.505726; ci = -0.00659119; ca = 4.09735; = -0.0270064, (6.79) 
which suggests that this stationary solution is unstable. 



When r > 0, the sign of C3 depends on whether one is taking the positive or negative sign in Eq. (6.63). For our 
initial conditions C3 is positive for /2+ and negative for /a- corresponding to the ± choice in Eq. (6.63). Combining 
these two equations we now find: 



5P =F ^'^SP = ji/) T ^"^H = 0, 



(6.80) 
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where 17^ = jciCsl > 0. This suggests that the fixed point associated with /2+ is unstable and stable under small 
perturbations. For the solution with r=0.05, for /i we find that 

ci = .034; C2 = 4.16915; n = .378701. (6.81) 

For /2+ 

ci = .034; C3 = 3.58; O = .351073. (6.82) 

For h- 

ci = .1251; C3 = -13.0305; 17 = 1.27676. (6.83) 
For the solution with r=0.005, we find for /i that 

ci = 0.00332219; C2 = 4.03982; n = 0.115849. (6.84) 

For /2+ 

ci = 0.00334441; = 3.95982; VL = 0.115079, (6.85) 

which suggests instability. For /2_ 

ci = 0.0406897; cg = -48.1771; n = 1.40011. (6.86) 

The small oscillation frequencies for the stable solutions /i and /2- are borne out by numerical simulations. 
However, for the predicted unstable case, at early times, < t < 650 the solutions exhibit a small oscillation with 
frequency Q, = 0.0448. After that time the solution switches to the separatrix solution with large oscillations in /3 and 
(j) but small oscillations in p and q. 

G. K = 1 

For comparison let us review what happened in the case n = \ which is quite different. At k = 1 the Lagrangian is 
given by: 

i = -2 r y^TT sech cos(0(i)) + 4^ {p{t)m - p{tf ~ m + 6 - + (6.87) 

From this we get the following four equations: 

|(/3p)=0, (6.88) 
/3 = -^y2g sechCsin(0(t)), (6.89) 
~ n~^.N 7rV2ff ?'tanhCsechCcos(0(t)) 

1 = -^p^^2 ' (6-90) 

< , x9 ~, . c ,9 x9 r fj(t)Tr'^^/2q tanhC sechCcosfoifi)) 
p(i)g(t) - p{t)' - 0(<) + <5 - fc2 + /?(t)2 ^ ^ ^ ^ (6.91) 

where 



From Eq. (6.88) we obtain 

l3{t)p{t) = constant — ai. (6.93) 
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Combining Eqs. (6.901 and (6.91) we find 



The stationary solutions have 



~2 r p{t)'K^^/2g tanhC sechC cos{(j){t)) 2 1/2 



+ k'- 



(j){t)=as\ p = Ps\ q = Vs; l3 = l3s; Cs = 

We have cf) = nir and we can restrict ourselves to 0, tt. Thus 

y/2g TT^ r tanh Cs sechCs 
Vs = 2ps T 



T^Ps 



(6.94) 



(6.95) 



(6.96) 



as =PsT 



_2 1" psTr^y/2g tanh Cs sechCs „2 



2/32 



(6.97) 



We are interested in small oscillations around the stationary solutions. We will choose f3 and cj) as the independent 
variables with p = ai//3, and C — 2fi^lt) ■ Letting l3 ~ j3s + Sj3, and cj) = a^t + 5(j) we obtain for small oscillations: 

(6.98) 



(5/3 = ^ — yj2gsec\iCs5(t> = TC2<5^, 



<50 = [213, 
= C3<5/3. 

Thus we get the equations 
When C = instead we have 



3 tanh Cs sechCs + sechCs (2sech^ Cs — 1] 

Ps 



Thus the stationary solution is 



(5/3 ± C2C35/3 = , (5(?!. ± C2C35<?!) = 0. 
P = 0;^ = 0, 

/? = -r7rV572sin(^), ^ = (5 - fc^ + /3^// 



And the small oscillation equations are 



/3s = 



5P = TrTT^/gj25^- 5^ = 2(3,513. 



5/3 



(6.99) 
(6.100) 

(6.101) 
(6.102) 

(6.103) 



VII. DAMPED AND FORCED NLSE (THEORY) 

In performing numerical simulations, one is interested in adding damping to the problem that we have studied 
earlier. The damped and forced NLSE is represented by 



d 

ot ax^ 



(7.1) 



where a is the dissipation coefficient. This equation can be derived by means of a generalization of the Euler-Lagrange 
equation 



d dC d dC dC _ dT 

dt dip^ dx dip* dip* dipl 



(7.2) 
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where the Lagrangian density reads 



-0* — re 



i(kx+e) 



and the dissipation function is given by 



Inserting the ansatz Eq. (5.1) into (7.4) we obtain 

F = -2a fl{Mt){x - qmiPiq -x)+pq- 



Integrating this expression over space we obtain 

F 



-2aCoy(M-0)- 



(7.3) 



(7.4) 



(7.5) 



(7.6) 



On the other hand L — J dxC is given by Eq.(5.17). Since F contains q and (j>, only the equations for f3 and p could 
be changed by the damping term. For k = 1, we know that the damping only affects the equation for f3 [S]. For 
K = 1/2 we have the same scenario, i.e. now the equation for /3 reads 



grC smjB) 
a"'^^ 6/3 sinh(C)' 



(7.7) 



The factor 2/3 comes from 2/(2/k — 1). 



VIII. NUMERICAL SIMULATIONS ON THE UNFORCED AND THE FORCED NLSE 

In this section we would like to accomplish three things. First, we would like to show that our numerical scheme 
as applied to exact soliton solutions with either r = or r ^ leads to known results (r — 0) and to new results that 
the exact solutions to the forced problem are metastable and at late times become a solitary wave whose parameters 
oscillate in time. Secondly, we want to compare the results of solving the four CC equations for the variational 
parameters /3, g, p, (p with a numerical determination of these quantities found by solving the FNLSE. We will find 
that if r is small compared with the amplitude of the solitary wave, that solution of the CC equations gives a good 
description of the wave function of the FNLSE. Finally, we would like to demonstrate that the regions of stability for 
the solitary waves of the FNLSE are well determined by studying the phase portrait as well as the p{v) curve for the 
approximate solution found by solving the CC equations. Most of the simulations will be restricted to k = 1/2 except 
for a discussion of blowup of solutions for k > 2. 



A. Numerical methodology 

Before displaying the results of our simulations, we would like to say a little bit about the numerical method we 
used and the boundary conditions, since we are constrained to perform the calculation in a box of length 2L. We 
have allowed the length to vary from 100 to 400. The number of points used on the spacial grid was 2L/Ax. The 
numerical simulations were performed using a 4th order Runge-Kutta method. + 1 points are used on the spatial 
grid n = 0, 1, . . . A'^. When studying the exact solutions we have used three different boundary conditions. For r = 
we use "hard wall" boundary conditions, where the wave function vanishes at the boundary: 

ip{±L,t) = 0. (8.1) 

For the case of r ^ and only when studying the exact solutions we have used mixed boundary conditions (see 



(4.27)). Otherwise, we have used periodic boundary conditions 

^P{-L,t)^i^{L,t), ^,{-L,t)^^j,{L,t). (8.2) 
In one of our simulations we have compared the use of periodic vs. mixed boundary conditions and found that the 



differences in the evolution of both /3{t) and ip{x,t) are hardly visible by the "eyeball method" (see Fig 12 ) The other 
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parameters related with the discretization of the system are increasing values of L, namely L = 50, L = 62.8, L = 100 
or L = 200. We have chosen as our grids in x and t, Aa; = 0.05 and At = 0.0001, (such that At < (l/2)(Aa;)2). 

We next need to determine the parameters q and /3 used in the CC equations. We determine g(t) from our numerical 
simulation by equating it to the value of x for which the density of the norm is maximum. 

To determine /3(t) for finite energy solitary waves, we assume that the variational parameterization of ip{x,t), 
namely Eq. (5.9) with A given by Eq. (5.10) is an adequate description of the wave function. If we do this and 
determine ijj{x = q{t)) then we have that 



Then we determine /3 as follows: 



\^P\\x = qit)) ^ A' 



1/k 



P^^[9n^W-{x^qm/{K + l). 



(8.3) 



(8.4) 



In the particular case when we are studying the time evolution of the exact solution for k — 1/2, i.e. Eq. (4.16), with 
conditions of Eq. (4.201, we need to subtract off the constant term to determine P{t). From Eq. (4.16) we can obtain 
P{t) from 
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One can also compute the momentum P{t) given by Eq. ( |2.12[ ). 

As an initial condition in most of our simulations (when we are not discussing the exact solution), 
an approximate solution given by the variational ansatz Eq. (5.9) with initial conditions for /3o, poj 9o 
In comparing with the CC equations, we will solve the four ordinary differential equations for (3, p, q, 
arbitrary k are given in Eq. (6.8) to Eq. (6.13). 



(8.5) 



we will use 
= and 4>o. 
(h which for 



B. Numerical simulations of PDE for r = 0, arbitrary k 



The initial conditions are those of the exact 1-soliton solution of the unforced NLSE given by (3.1 ). For r = and 
(5 = 0, the stability of the NLSE has been well studied as a function of k. For k < 2 the solutions are known to be 
stable, and for k > 2 the solutions are unstable. For k = 2 there is a critical mass above which the solutions are 
unstable. The nature of the solution when it is unstable has been studied in variational approximations [19^ as well as 
using various numerical algorithms |29l 132] . Here we want to show that our code reproduces these well known facts. 
We are also interested in the effect of S in our simulations and we find that the critical value of k is independent of 
S. This agrees with our arguments earlier on the effect of scale transformations on the stability of the exact solutions 
for r = 0. That is, stability does not depend on S. 



C. Simulations at r = 0, different values of k 

In this section we study the numerical stability of the exact solutions for r = for different values of k (k e [0.25, 3]). 



First we set r = 0, (5 = 0, 5 = 2 in NLSE and we start from the exact soliton solution (3.1 ) with /3 = 1, p = 0.01 and 
(j>o — 0. In this simulation the boundary conditions chosen were ^{zLL,t) = 0. We notice for k < 2 as shown in Fig. 
^ the solitary wave moves to the right and maintains its shape. 

Indeed, in the simulations the solitons are not stable for k > 2 (see Fig. [?]), the amplitude of the unstable soliton 
grows and the soliton becomes narrow. Also the soliton moves only slowly to the right while blowing up at a finite 
time. In Fig. [8] upper panel we show the real and imaginary parts of the field for k — 1.5,2,3 for the final time 
of the simulations. The result of adding a term proportional to S does not alter the behavior of the solitary waves. 
Choosing 5 = —1 we obtain similar results as for S — 0, i.e. the soliton is unstable for k > 2. In Fig. [8] lower panel 
we show the evolution of /3 for short times and k ~ 1.5; 2; 3. This shows once we reach the metastable solution for 
K = 2, we see P{t) 00 in a finite amount of time. This agrees with our analysis based on Derrick's theorem (see 



Eq.(3.12 )), and with the discussion of the critical mass needed for blowup for k = 2 in [TS]. For S = —I we have also 
investigated the constant phase solutions that fulfill the condition (3.5). We have studied numerically the case r = 0, 
S = —1, g = 2 with initial conditions p — 0.01, t/fQ = and (3 = ^—k'^{S + p^) and varied k £ [0.25, 3]. We found that 
again the soliton becomes unstable for k > 2 in accord with our result that instability as a function of k does not 
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depend on S. Our numerical experiments for the case r = show that our codes reproduce well known results for the 
stability of the solutions. In a future paper we will compare the numerical solutions in the unstable regime with the 
predictions of the CC method. Here we will use the exact form of the solution rather than the post-Gaussian trial 
functions we used earlier [19 as well as include an additional variational parameter in the phase that is canonically 
conjugate to the width parameter /3. To estimate the finite size effect on the definition of /3 we notice in Fig. [9j that 
the values of /3{t) in the simulation of the NLSE deviate from the exact constant value of 0.5 by approximately 0.25%. 
By increasing L one could reduce this error. 



D. On the exact solutions for r 7^ and k = 1/2 



Earlier we showed that for r 7^ and k ~ 1/2 there are exact solutions of the form (see Eq. (4.24)) 

2 



fiy) = -sign(r) [a + bsech^(3{y, t)] , 

(8.6) 



with a, /3 and b given by Eq. (4.20). In the numerical simulation shown in Fig. |9]and Fig. 10 the parameters used in 



the simulation are k = 1/2,5 — —l,g — 2,k = 0.1, r = —0.075. For that case: 

f{y) = 0.486113sech^(0.402539y) + 0.0904622. (8.7) 

The unstable stationary variational solution corresponding to this exact solution is 

/2+(y) = 0.745535sech2(0.498509y). (8.8) 

These solutions are shown in Fig. [2] As we can see from the numerical solution of the FNLSE shown in Fig. [9] the 
exact solitary wave oscillates in amplitude and width which shows that the original exact solution was unstable. In 
spite of this the solitary wave neither dissipates nor blows up. The width parameter oscillates from its initial value 
of around 0.4 to an upper value of around 0.8 with a period T of around T = 20. If one reduces the driving terms 
to fc = 0.01 and r = —0.005, the solitary wave again oscillates in amplitude and width but the amplitude of the 
oscillation is only 10% of the total amplitude and the oscillation frequency is reduced from the previous case as is 
shown in Fig. \TT\ 



In Fig. 12 we compare the behavior of /?(<), and the late form of the wave functions using two types of boundary 
conditions (BCs): periodic BCs {tp{—L,t) = tp{L,t), %px{—L,t) = ij}x{L,t)) shown as dashed lines, and "mixed" BCs 
which are shown as solid lines. In this simulation the two boundary conditions lead to almost identical results. Most of 
the simulations (except the ones starting from the exact solution) were performed using periodic boundary conditions. 
Next we look at the case where r > where for amplitudes A{t) > there are no exact solutions. However, the 
stationary solutions of the CC equations that are stable to linear perturbations do lead to solitary waves that have 
widths whose oscillations have only very small amplitude. To show that we start with the stationary solutions with 
different positive r but other parameters (k, 5, g, k) the same as in Fig. [9j We show this in Fig. 13 where we consider 



solitary waves moving to the left initially with constant velocity v = —2k. We see that as we increase r the amplitude 
of the solitary wave increases and in Fig. [14] we see that the average value of f3 also increases with r. We notice in 
Fig. [14] that the oscillations in f3 get more chaotic as we increase the value of r. 



E. Comparison of numerical simulations of the FNLSE with the solution of the equations for the collective 

coordinates 



In this section we would like to show that the solution to the CC equations gives quantitatively good results for 
the collective variables p{t), f3{t),q{t) and (j){t) when these quantities are calculated directly from the solution of the 
FNLSE. We also want to show that the criterion for stability of the CC equations, namely a study of the phase 
portrait or the p{v) curve, gives an accurate measure of what initial values of collective variables lead to stable vs. 
unstable solitary wave solutions. First we discuss the linearly stable stationary solution Eq. (6.70) that is shown in 
Fig. [5 In the CC equation (3 remains at the fixed value /3(0) = /3s = 0.51773. Using this as an initial condition in the 
FNLSE solved using periodic boundary conditions, one finds that the solitary wave has a mean value of /3 only 2 % 



different from the solution, Fig. 15 For < t < 30, phonons (short for linear excitations) are radiated by the soliton 
which travel faster than the soliton. Coming from a boundary these phonons reappear and collide with the soliton 
producing the increased oscillations seen in Fig. [l5| For the system of length 2L — 125.6, these increased oscillations 
begin at tc « 125. When the system size is doubled, both the soliton and the phonons have to cover doubled distances 
before the collisions begin at tc ~ 250. 



28 



If we add a little bit of damping and start from a value of /3q which dissipates to a stationary point then all the 
collective variables are well approximated by the solution of the CC equation as shown in the middle panel of Fig. 
[T6l Here we have chosen g = 2, k = 0.1, r — 0.05, S = —1, 9 = 0, a = 0.05 . For these values, the stationary solution 
is approximately given by Eq. (6.69) which has a value of /3s — 0.50512. In the middle panel of Fig. 16 we use the 
initial condition — 1. We find both the CC equations and the solution of the FNLSE relax to the approximate 
stationary solution of Eq. ( |6.69 ) . 

In Fig. [it] we turn off the damping and see how both the simulation and CC equations evolve. Solving the CC 
equations with the initial condition /3q = 1, one finds from both the orbit of this initial condition in the phase portrait 
(see left bottom panel) and the p{v) curve that the solitary wave should be stable. In the middle panel we see that 
both the exact solution and the solution to the CC equation for /? oscillate with 1% oscillations around either the 
initial value /3o = 1 (for the numerical solution) and a slightly shifted value (0.994) for the CC equation. The actual 
solution has another oscillation frequency for the height at maximum. 



In Fig. 18 we study the time evolution of a stationary solution that is known to be unstable. We use the of solution 
Eq.(6.71). Both the CC equations and the simulation of the FNLSE show that this solitary wave develops into a 
solution that is represented by a separatrix in the phase portrait. The negative slope of the p{v) curve also predicts 
instability. 

The soliton stability depends strongly on f3o. For the parameters and initial conditions of Fig. 19 both stability 
criteria predict the following pattern: instability for /3o < 0.484511, stability for 0.484511 < 13q < 0.5458, instability 
for 0.5458 < f3o < 0.65 and stability for /3q > 0.66. In Fig. 19 we present two examples: a stable soliton (/3o — 0.53) 
and an unstable one (/?o — 0.65), both confirmed by our simulations. We note, however, that the above boundaries 
between stable and unstable regions do not always fully agree with the simulations: the errors vary between 1.7% and 
12%. 

In Fig. [20] we show a particular case where the time evolution of an initial condition which was close to an unstable 
stationary solution of the CC equation exhibits intermittency in both the inverse width (and thus the amplitude A 
also) as well as the energy. We have also observed intermittency in the solutions of the CC equations for k = 1. 
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FIG. 6: r = 0. Soliton moving to the right at t* = 166.6, 333.3, 500, 666.6, 833.3, 1000. Left and right panel: k = 0.5 and 
K = 1.5, respectively. Parameters: 5 = Q, g = 2 with initial conditions /3 = 1, p = 0.01 and (f)o = 0. 



IX. CONCLUSIONS 



In this paper we have studied analytically, numerically and in a variational approximation the FNLSE with arbitrary 
nonlinearity parameter k. We studied in detail a variational approximation based on the solutions to the unforced 
problem and have studied the CC equations coming from the Euler-Lagrange equations. We determined the stationary 
solutions of the CC equations and studied their linear stability. We found that for small forcing parameter, the CC 
equations give quantitative agreement with directly solving the FNLSE. Also the domains of stability of the initial 
conditions for the FNLSE were quite remarkably close to those found by studying the orbits in the phase portrait 
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FIG. 7: r = 0. Solitary wave moving to the right (blowup when k > 2). Left panel: k = 2.0 (at t* = 
166.6; 333.3, 500, 666.6, 833.3, 1000), respectively. Right panel: k = 2.25 (at t* = 5.8; 11.6, 17.5, 23.3, 29.2, 35.0). Parameters: 
6 = 0, g = 2 with initial conditions /3 = 1, p = 0.01 and (f>o = 0. 
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FIG. 8: r = 0. Real (solid line) and imaginary (dashed line) parts of the fields for k = 1.5 for t* = 1000 (left upper panel); 
K = 2.0 for t* = 120 (right upper panel); n = 3.0 for t* = 25 (left lower panel). Evolution of /3 from the direct simulations of 
NLSE. K — 1.5 (solid line); k = 2 (dashed line); /t = 3 (dotted line) (right lower panel). Parameters: 5 = 0, g = 2 with initial 
conditions P = 1, p = 0.01 and (po = 0. 



for the CC equations as well as the stability curves p{v). We also found that the linearly stable stationary solutions 
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FIG. 9: Test of the exact solution Eq. (4.161 (with the parameters a, /3 and b determined by Eq. (4.19|) of the perturbed NLSE 
for K = 1/2, S = —1, g = 2, k — 0.1. Here we display the inverse width parameter /3(t) of the soliton. Left panel; simulation of 
the unperturbed NLSE with r = (a = 0) corresponding to the constant phase solution of Eq. (3.6 1, P'^ = —5/A, Right panel: 
simulation of the perturbed NLSE with r — —0.075. 
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FIG. 10: Soliton energy computed by subtracting the background energy density from the total energy density. Left panel: 
r = (a = 0). Right panel: r = —0.075. Parameters: are the same as in Fig. [9] 



of the CC equations were quite close to the stationary solutions of the FNLSE, and when these stationary solutions 
were used as initial conditions for the FNLSE the oscillations of the width parameter were of small amplitude. These 
simulations at k = 1/2 reinforce our belief that our two dynamical criteria for understanding the stability of the 
solitary wave, namely the stability curve p(v) as well as studying orbits in the phase portrait are accurate indications 
of the stability of the solitary waves as obtained by numerical simulations of the FNLSE. We intend to extend our 
study of the FNLSE equation to the regime k > 2 to see how forcing affects the blowup of solitary wave solutions 
in that regime. This will be done in the variational approximation by adding a variational parameter canonically 
conjugate to the width parameter f3 similar to the approach of (TH] [T^. Our results are likely to shed light on the 
behavior of a number of physical systems varying from optical fibers [1] , Bragg gratings [2] , BECs [3l |4] , nonlinear 
optics and photonic crystals [7]. 
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FIG. 11: Time evolution of the exact solution of the FNLSE as in Figs. |9] and [To] but using smaller value of the driving terms: 
k = 0.01 and r = -0.005. 
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FIG. 12: Comparison of the two boundary conditions: periodic (dashed line) vs mixed (solid line) r — —0.075. Left panel: 
I3{t), right panel: soliton profile at t* — 200. The parameters are the same as used in Fig. [9] and Fig. 10 
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FIG. 13: Here we compare the time evolution of a solitary wave traveling to the left initially with constant velocity — 2fc 
starting from an approximate stable stationary solution for three different values of positive r. The times are f = 33.3, 66.6, 
100, 133.3, 166.6, 200. We use mixed boundary conditions. Left panel: r = 0.01, middle panel: r = 0.05, right panel: r = 0.075. 
Parameters: k = 1/2, 5 = —1, g = 2, k = 0.1. Initial condition (|4.16| with (|4.19K 
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FIG. 14: Time evolution of the inverse width parameter P{t). Here we plot /3(t) for the same initial conditions as in Fig. 13 
Again we have left panel r = 0.01, middle panel r = 0.05 and right panel r = 0.075. 
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FIG. 15: Test of stable stationary solution of the CC equations (horizontal line) by a simulation of the NLSE with periodic 
boundary conditions (black and red solid lines for system sizes 125.6 and 251.2, respectively, the latter line is shifted upwards 
by 0.05). The influence of linear excitations, which are radiated at early times, on the soliton oscillations is discussed in the 
text. The parameters used are k = 1/2, r = 0.05, k — —0.1, 9 — 0, g — 2, S — —1, a = 0, with initial conditions po — —k, 
qo = 0, (j>o = and fio — jSs = 0.51773. For the CC equations we choose po = — + lO^'' to avoid numerical singularities. The 
approximate solitary wave is described by Eq. (6.701, and is shown in Fig. [5] 
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FIG. 16: Comparison of simulations of the collective coordinate (CC) approximation and the numerical solution of the PDE 
where we have considered the case k = 1/2 and included damping when appropriate. Left upper panel; Soliton moving to the 
left for t* — 250,500. Right upper panel: Real (solid line) and imaginary (dashed line) parts of tl>ix,t) for t* = 500. Real 
(red solid line) and imaginary (blue dashed line) parts of the exact solution for the undamped NLSE for t = 500. Left and 
right middle panels: comparison of the time evolution of /3 and q computed from the simulations of the PDE (solid line) and 
numerical solutions of the CC equations (dashed line). Lower panels (different scales are shown): for t* = 500 soliton profile 
from simulations (red solid line) and from the exact solution (4.161 and (4.191 (red dashed line). Notice that the exact solution 
is obtained for a — 0. Parameters: g = 2, k = 0.1, r — 0.05, 5 = —1, 9 = 0, a = 0.05, with initial conditions /3o = 1, po = 0, 
go = and (f>o = —1.69 + tv/2 ~ —0.119. Notice that for this set of parameters the condition |r| < fc'^/(4gr) (see below Eq. 
(4.19 1 ) is satisfied. After a transient time, since the numerical solution approaches the stationary solution, the p{v) curve is 
just a point (-0.1, -0.2). The phase portrait is also represented by a point (0.04491169724, -0.02819002364). 
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FIG. 17: Results from simulating the time evolution of the solitary wave using both the PDE as well as the CC equations for 
n = 1/2 with no damping. We use the same parameters as in Fig 16 but a = Q. Upper panels: soliton profiles (different 
scales) for t — 250 (black) and t = 500 (red). Middle panels: time evolution of /3 computed from the simulations of PDE (left) 
and computed from the numerical solutions of CC equations (right). Lower left panel: elliptic orbit in the phase portrait with 
positive sense of rotation which predicts stability; lower right panel: positive slope of the p{v) curve predicts stability. 
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FIG. 18: Unstable stationary solution for 
solution has po = ~k (simula tions) and po 
(see Eq, 



Vfc'* - 8rg/3)/8) 



1/2 



0.05, K = 1/2, 5 = 2, fc = -0.1, 6* = 0, (5 = -1 and a = 0. The stationary 
— k + 10~"' (numerical solutions of the CCs), cj>o = n and /3o = /9s = {{k'^ + 



6.711. Here we find for P{t) that both the simulation (solid line) and the solution to the CC 
equations (dotted line) show switching to a new solution. From the p{v) curve the negative slope of the curve predicts this 
instability. In addition, the orbit of the phase portrait is a separatrix, which also predicts the instability seen in /3. 
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FIG. 19: Simulations for the NLSE with r = 0.05, k = —0.1, periodic BC, with initial conditions /3o — 0.53 (black curves) and 
0.65 (red). The other parameters and initial conditions are k = 1/2, g = 2, 5 = —1, 9 = 0, a = 0, (l>o = 0, po ~ —k + 0.001 
and go = 0. Upper left panel: P{t) compared with numerical solutions of the CC equations (dotted lines). Upper right panel: 
orbits of the phase portrait. Lower panels: p{v) curves. 
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FIG. 20: Numerical simulation of the NLSE for r — 0.05 using periodic boundary conditions and k — 1/2. /3(f) and E{t) show 
intermittency. Upper panels: Soliton moving to the right for t* — 1000 (black) and t* = 2000 (red) and the evolution of the 
energy of the system. Lower panels: time evolution of /3 and q computed from the simulations of PDE. Other parameters of the 
simulations: g = 2, k — —0.01, S — —1, — 0, a = 0, with ini tial co nditions pp = Ps — —k, go = 0, (jio = ti" and Po — Ps = 0.4983. 
This corresponds to the unstable stationary solution of Eq. (6.731, namely /2+(y) = 0.745042sech^(0.498345j/). 



